Derivation of the residence time for kinetic Monte Carlo simulations 
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The kinetic Monte Carlo method is a standard approach for simulating physical systems whose 
dynamics are stochastic or that evolve in a probabilistic manner. Here we show how to calculate 
the system time for such simulations. 
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I. INTRODUCTION 

In a kinetic Monte Carlo (kMC) simulation, a model 
system is moved from state to state according to prob- 
abilistic rules. A nice example is adatom diffusion on a 
surface: in this case the "system" is the adatom plus 
surface, and a "state" has the adatom at a particu- 
lar surface site. The probabilistic rules for moving the 
adatom from one site to another (i.e., the jump rates) are 
obtained from molecular dynamics (MD) calculations, 
which give the energy barriers over which the moves oc- 
cur. Subsequently, a separate kMC simulation moves the 
adatom from site to site over the surface, and thereby 
enables calculation of the adatom diffusion coefficient, 
D = ^a;^) /(4i), from distance x traveled by the adatom 
over the time interval t. 

This example is typical in that the kMC simulation 
gives a macro- or mesoscopic quantity of interest (in this 
case the adatom diffusion coefficient), using information 
obtained from the microscopic or atomic scale. In fact, 
molecular dynamics calculations alone are impractical for 
simulating adatom diffusion over a surface: the adatom 
at room temperature attempts a jump every 10"^'^ sec- 
onds, while it is successful once every 10^^ seconds (so in 
MD simulations it is the numerous unsuccessful attempts 
at many sites, rather than the rare successful one, that 
enables the energy landscape to be mapped and jump 
rates determined). 

The great virtue of the kMC method is that it en- 
ables simulation of physical processes that involve very 
disparate time scales. It is then essential that the actual 
time in the physical system be correctly reproduced by 
the residence times in the states of the model system, 
or equivalently, by the time between events in the model 
system. 

While equivalent statements, the different semantics 
suggests the ability of the kMC method to obtain prop- 
erties of an equilibrium system as well as to model the 
dynamic evolution of a non-equilibrium system. In the 
case of equilibrium systems, each state may be visited 
many times, so that the residence times in the states 
are proportional to the equilibrium populations of the 



physical states (for example, the adatom on a finite sur- 
face may be found at sites in proportion to the residence 
times at the corresponding states in the model system). 
More generally the kMC method is used to model the 
evolution of a non-equilibrium system (that is, where a 
state once left can never be returned to) , so this case will 
motivate the derivation of the residence time presented 
below. However, the mathematics apply to equilibrium 
systems as well, as will be made apparent at the end of 
the derivation. 



II. DERIVATION OF THE RESIDENCE TIME 

Consider a collection {j} of elements, where each el- 
ement j has an expected lifetime Tj. The kMC method 
follows from the assumption that the values {r,} are 
time-independent (i.e., that the transition rates {t^^} 
from the current state to accessible states are time- 
independent). Then the probability pj{t)dt that the par- 
ticular element j will fail during the infinitesimal time 
interval [t, t -f dt] is given by the probability that element 
i will not fail prior to time t, multiplied by the proba- 
bility T^^dt that it will fail during the subsequent time 
interval dt; that is, 



Pj{t)dt= 1-Jpj(t')dt' 



This equation simplifies to 

t 



1 



Pj{t')dt' 



dt 



(1) 



(2) 



Taking the derivative of each side of Eq. ^ with respect 
to t and integrating produces the exponential probability 
density function (PDF) 



1 



Pj(t) = — exp 



(3) 
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Specifically, Pj{t) is the distribution of failure times for 
element j; the average value of the failure time is (tj) = 

J^tp-i(t)dt ^ Tj. ^ 

To perform a simulation, we need to randomly select a 
time t from this distribution. The formula for converting 
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a random number x taken from the uniform probability 
distribution P{x) = 1 (such x values are produced by the 
standard random number generators) to the correspond- 
ing t value is derived as follows. The probabilities p{t)dt 
and P{x)dx must be equal, so p{t)dt — dx. Thus 



x{t) = j p{t')dt' = 1 - exp 





t 



(4) 



Inverting this expression then gives the desired relation 
between x and t, 



t = Tj[- ln(l - x)] 



(5) 



with X randomly chosen from the interval [0, 1). 

In this way we get the set {tj} of element failure times. 
If the elements {j} are completely independent of one 
another, they fail according to the time ordering of the 
set {tj]. 

If, however, the elements are not independent, the re- 
maining values {tj} are altered with each failure. Con- 
sider that the smallest member of the set {tj] is Ti, so 
the first element to fail does so at time Ti. Which is 
the next element to fail, and when docs that occur? If 
element j is still functioning at time Ti, the PDF 
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Pjit ~ Ti) = — exp 



(6) 



where t' is the new lifetime of element j for time t > T\. 
Equation ([6]) gives the distribution of times i > Ti at 
which element j fails, or equivalently, the distribution of 
time intervals At = t — Ti at the end of which element 
j fails. As above, a value At may be randomly chosen 
from this distribution by use of the relation 

At = Tj[-ln(l - x)] (7) 

where x G [0, 1). Thus the next element to fail is that 
producing the smallest member of the set {Atj}. 

An alternative approach to simulating the system evo- 
lution is to randomly choose the next element to fail 

according to the set of probabilities ^t^^ ('^/^)}' 

where element k is one of the set {j}, t^T^/X), i'^i^) 
the probability that element k is the next to fail, and the 
sum is over all N currently surviving elements j (remem- 
ber that the values {tj} may change after every failure). 
This failure occurs at the end of the time increment At, 
that may be taken to be the average value of the smallest 
member of the set {Atj} (were the set to be calculated 
inmmrerable times). 



TV 
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Atk 
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k=l 



1 Efe[-Ml-2:fe)] 



N 



E 



ln(l - Xk)) 



(8) 



A different expression for At may be found by noting 
that the probability p{At)dt that the next element to fail 
will do so during the infinitesimal time interval [At, At + 
dt] is 



N 



N 



p{At)dt^J2lpj{At)dt Y[ [1-pkiAt)] 



(9) 



fe(#j)=i 



where the content of the curly brackets is the probabil- 
ity that element j will fail during the infinitesimal time 
interval [At, At -I- dt] , multiplied by the probability that 
no other element will fail during At. Then 



p{At) = 51 ^ - 

.7 = 1 '^^ 



At 

Tn 



N 



n '^^p 



At' 



N 

j = l J 



■ exp 



N 



gr-M exp ^-AtE^-L 



'3 = 1 3 



(10) 



As above, a value At may be randomly chosen from this 
distribution by use of the relation 



At 



-ln(l -x) 



(11) 



where x £ [0, 1). The average value is given by Eq. ([5]), 
as expected. 

The evolution of the system of elements is thus accom- 
plished by selecting an element to fail, and incrementing 
time accordingly. With each failure the system enters a 
new state; obviously it cannot return to any old states. 
For this reason it has been convenient to use the set of 
"lifetimes" {tj} corresponding to the surviving elements 
J, rather than the set of transition rates {ki^j{ for tran- 
sitions from the current state i to the accessible states j 
(the system will enter state j if element j is the next to 
fail). The connection between the sets is made by noting 
that Ej''"j~^ ~ 'Ylij^i^j^ where the system is currently in 
state i in either case. 

By making this replacement in the denominators of 
Eqs. ([5]) and pTjl . those equations for At may be used 
for equilibrium systems (where the set of states and tran- 
sition rates between states don't change) as well. More 
generally, Tj in the equations above may be replaced by 
^hen it is understood that the system is currently 
in state i. 

For computational efficiency, it is most typical for a 
simulation to select the transition from current state i ac- 
cording to the set of probabilities |fci-.j'/Ej^j^j | ^'^d 
calculate the transition time interval At using either Eq. 
([8]) or pT|) . rather than to calculate the set {At^} from 
the relation Atj — k~2,j[— ln(l — x)] and choose the des- 
tination state j' from the smallest of those values. When 
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the distribution of simulation completion times is desired 
(for example, the time at which all elements have failed), 
it is necessary to perform a large number of nominally 
identical simulations, using Eq. pT|) rather than Eq. ^ 
for every transition time interval At. (Note that this 
distribution will converge to a Gaussian distribution, ac- 
cording to the central limit theorem.) 

More generally, evolution of a dynamic system is not 
described by a set {fci— j } of time-independent rate con- 
stants. In this case there is no set of probabilities 

|fci_>j'/^^A;i^j I from which to select the destination 

state j' . However, a set {tj} of transition times can still 
be calculated from a PDF, giving the sequence of events 
from an initial time. For example, consider that the dis- 
tribution of failure times for element j is given by the 
two-parameter (7 and r) WeibuU distribution 



p{t) = 7 



■ exp 



71 



(12) 



rather than by Eq. ^ . For 7 > 1 this PDF decays faster 
than exponentially and so might be appropriate for an 
element that "ages", or accumulates damage over time. 
(Indeed, the probability TT(t)dt that the clement, having 
survived to time t, will immediately fail is j{t'^^^ /T^)dt, 
meaning that the probability of immediate failure in- 
creases with time.) The first and second moments of 
this distribution are, respectively. 



and 



7 V7 
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(13) 



(14) 



where F is the Gamma function. The counterpart to Eq. 
O is 



x{t) = J p{t')dt' = 1 - exp 



so the failure time for element j is 

i,=r,[-ln(l-x)]i/^^ 



(15) 



(16) 



where x G [0, 1). 

If the parameter values {7^ , Tj } change with each ele- 
ment failure (reflecting, say, more rapid aging), the PDFs 
{pj{t)} given by Eq. (fT^ change accordingly. To calcu- 
late a new set {tj} of transition times following an ele- 
ment failure at time T, we note that x and t must be 
related by 



dx' 



p{t')dt' 



(17) 



where the lower limit x(T) = 1 — exp[— (T/r')''' ] is ob- 
tained from Eq. (|15p . and 7' and r' are the parameter 
values for t > T. Then the counterpart to Eq. (jl]) is 



x(t) = x(T) + exp 



1 — exp 



y \ 7 

r' 







— exp 









(18) 



(19) 



so the failure time for element j is 

i,=r;[-ln(l-x)]i/^i 
where x G [x{T), 1). 

III. DISCUSSION 



In physics and materials research, the kMG method is 
most often used for simulation of non-deterministic, ther- 
mally activated processes, as these give time-independent 
rate constants (time-independent probabilities) for tran- 
sitions between states. A similar, but more complex, 
example to the adatom diffusion calculation mentioned 
in the Introduction is diffusion of vacancies and intersti- 
tial atoms in a grain boundary [l| , where the transitions 
between states may occur by very complex mechanisms 
(for example, by concerted motion of atoms). However 
the rates k are of the simple form 



ki^j = V exp 



k^ 



(20) 



where E^, is the energy barrier between the initial state 
i and the final state j (which the system may overcome 
with thermal energy supplied by local temperature fluc- 
tuations), /cb is the Boltzmann constant, T is the average 
temperature of the system, and v is the frequency with 
which the system attempts to make the transition. The 
latter quantity gives the time scale for the microscopic 
process, that is needed for the kMG simulation. 

Besides its natural application to atomic-scale pro- 
cesses, the kMG method has been used to simulate the 
time-dependent damage and failure of material under 
stress. For example, Gurtin et al. [3] considered a spring 
network model where the rate r of failure at site j is 



r,{t)^Aa,{tr 



(21) 



x{T) 



where (jj (t) is the local stress at time t, and the exponent 
Tj > 1 ensures a nonlinear relationship between damage 
rate and stress as is the case for creep. The time between 
failures was calculated by Eq. ([8]), after which all the lo- 
cal stresses {aj} and failure rates {rj} were recalculated. 
In contrast, Harlow et al. ^ considered a polycrystal 
under stress, and randomly chose, from a probability dis- 
tribution that included the local stress as a parameter, a 
time-to- failure tj for each grain boundary facet j. Then 
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time was advanced by the amount of the smallest mem- 
ber of the set {tj}, the load supported by that failed 
facet was redistributed to adjacent facets, and a new set 
{tj} was obtained from new distributions. In a similar 
way Andersen and Sornette considered a fiber rupture 
model where failure of each fiber affected the failure-time 
probability distributions for the remaining fibers. 

Effective use of the kMC method obviously requires 
that all important states of the system, and all important 
transitions between them, be identified. How to ensure 
this, and other technical and computational (i.e., practi- 
cal) issues are addressed in a review by Voter Q . Another 
problem is presented by "basins" of states, which are sub- 
sets of mutually accessible states from which escape is a 
very rare event. This of course defeats the purpose of 
the kMC method. However Van Siclen 6] has recently 
shown how to calculate the residence time for a system 
trapped in a basin, under the assumption that the sys- 



tem has equilibrated in the basin (which is to say, there 
is no correlation between the entry and exit states). 

To conclude, the kinetic Monte Carlo method is a pow- 
erful technique for simulating the dynamics or evolution 
of non-deterministic systems. Thus it should have ap- 
plications beyond the physical sciences, for example to 
biology, ecology, risk assessment, and finance. 
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